library(ggplot2)

exp1 <- read.table(file="exp_noed_nohet_maxev2.txt", header=TRUE)
exp2 <- read.table(file="exp_noed_nohet_maxev3.txt", header=TRUE)
exp3 <- read.table(file="exp_noed_nohet_maxev4.txt", header=TRUE)
exp4 <- read.table(file="exp_noed_nohet_maxev5.txt", header=TRUE)
exp5 <- read.table(file="exp_noed_nohet_maxev6.txt", header=TRUE)
exp6 <- read.table(file="exp_noed_nohet_maxev7.txt", header=TRUE)
exp7 <- read.table(file="exp_noed_nohet_maxev8.txt", header=TRUE)
exp8 <- read.table(file="exp_noed_nohet_maxev9.txt", header=TRUE)
exp9 <- read.table(file="exp_noed_nohet_maxev10.txt", header=TRUE)

msqer11 <- mean((exp1$betafrailgap-(-1))^2, na.rm=TRUE)
msqer21 <- mean((exp2$betafrailgap-(-1))^2, na.rm=TRUE)
msqer31 <- mean((exp3$betafrailgap-(-1))^2, na.rm=TRUE)
msqer41 <- mean((exp4$betafrailgap-(-1))^2, na.rm=TRUE)
msqer51 <- mean((exp5$betafrailgap-(-1))^2, na.rm=TRUE)
msqer61 <- mean((exp6$betafrailgap-(-1))^2, na.rm=TRUE)
msqer71 <- mean((exp7$betafrailgap-(-1))^2, na.rm=TRUE)
msqer81 <- mean((exp8$betafrailgap-(-1))^2, na.rm=TRUE)
msqer91 <- mean((exp9$betafrailgap-(-1))^2, na.rm=TRUE)


msqer12 <- mean((exp1$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer22 <- mean((exp2$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer32 <- mean((exp3$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer42 <- mean((exp4$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer52 <- mean((exp5$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer62 <- mean((exp6$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer72 <- mean((exp7$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer82 <- mean((exp8$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer92 <- mean((exp9$betafrailelapse-(-1))^2, na.rm=TRUE)


msqer13 <- mean((exp1$betaag-(-1))^2, na.rm=TRUE)
msqer23 <- mean((exp2$betaag-(-1))^2, na.rm=TRUE)
msqer33 <- mean((exp3$betaag-(-1))^2, na.rm=TRUE)
msqer43 <- mean((exp4$betaag-(-1))^2, na.rm=TRUE)
msqer53 <- mean((exp5$betaag-(-1))^2, na.rm=TRUE)
msqer63 <- mean((exp6$betaag-(-1))^2, na.rm=TRUE)
msqer73 <- mean((exp7$betaag-(-1))^2, na.rm=TRUE)
msqer83 <- mean((exp8$betaag-(-1))^2, na.rm=TRUE)
msqer93 <- mean((exp9$betaag-(-1))^2, na.rm=TRUE)


msqer14 <- mean((exp1$betacondgap-(-1))^2, na.rm=TRUE)
msqer24 <- mean((exp2$betacondgap-(-1))^2, na.rm=TRUE)
msqer34 <- mean((exp3$betacondgap-(-1))^2, na.rm=TRUE)
msqer44 <- mean((exp4$betacondgap-(-1))^2, na.rm=TRUE)
msqer54 <- mean((exp5$betacondgap-(-1))^2, na.rm=TRUE)
msqer64 <- mean((exp6$betacondgap-(-1))^2, na.rm=TRUE)
msqer74 <- mean((exp7$betacondgap-(-1))^2, na.rm=TRUE)
msqer84 <- mean((exp8$betacondgap-(-1))^2, na.rm=TRUE)
msqer94 <- mean((exp9$betacondgap-(-1))^2, na.rm=TRUE)


msqer15 <- mean((exp1$betacondelapse-(-1))^2, na.rm=TRUE)
msqer25 <- mean((exp2$betacondelapse-(-1))^2, na.rm=TRUE)
msqer35 <- mean((exp3$betacondelapse-(-1))^2, na.rm=TRUE)
msqer45 <- mean((exp4$betacondelapse-(-1))^2, na.rm=TRUE)
msqer55 <- mean((exp5$betacondelapse-(-1))^2, na.rm=TRUE)
msqer65 <- mean((exp6$betacondelapse-(-1))^2, na.rm=TRUE)
msqer75 <- mean((exp7$betacondelapse-(-1))^2, na.rm=TRUE)
msqer85 <- mean((exp8$betacondelapse-(-1))^2, na.rm=TRUE)
msqer95 <- mean((exp9$betacondelapse-(-1))^2, na.rm=TRUE)


msqcondfrail <- c(msqer11, msqer21, msqer31, msqer41,msqer51,msqer61,msqer71, msqer81, msqer91)
msqag <-  c(msqer13, msqer23, msqer33, msqer43,msqer53,msqer63,msqer73, msqer83, msqer93)
msqfrailelapse <- c(msqer12, msqer22, msqer32, msqer42,msqer52,msqer62,msqer72, msqer82, msqer92)
msqcondgap <- c(msqer14, msqer24, msqer34, msqer44,msqer54,msqer64,msqer74, msqer84, msqer94)
msqcondelapse <- c(msqer15, msqer25, msqer35, msqer45,msqer55,msqer65,msqer75, msqer85, msqer95)
      
sample <- c(2,3,4,5,6,7,8,9,10)

yvalues <- c(msqcondfrail,msqcondgap,msqcondelapse,msqag,msqfrailelapse)
xvalues <- c(rep(sample,5))
names <- c(rep("Conditional Frailty",9), rep("Conditional, Gap",9), rep("Conditional, Elapse",9), rep("Andersen-Gill",9), rep("Frailty, Elapse",9))
noednohet <- data.frame(msqerror=yvalues,sample=xvalues,Model=names)

exp1 <- read.table(file="exp_ed_nohet_maxev2.txt", header=TRUE)
exp2 <- read.table(file="exp_ed_nohet_maxev3.txt", header=TRUE)
exp3 <- read.table(file="exp_ed_nohet_maxev4.txt", header=TRUE)
exp4 <- read.table(file="exp_ed_nohet_maxev5.txt", header=TRUE)
exp5 <- read.table(file="exp_ed_nohet_maxev6.txt", header=TRUE)
exp6 <- read.table(file="exp_ed_nohet_maxev7.txt", header=TRUE)
exp7 <- read.table(file="exp_ed_nohet_maxev8.txt", header=TRUE)
exp8 <- read.table(file="exp_ed_nohet_maxev9.txt", header=TRUE)
exp9 <- read.table(file="exp_ed_nohet_maxev10.txt", header=TRUE)

msqer11 <- mean((exp1$betafrailgap-(-1))^2, na.rm=TRUE)
msqer21 <- mean((exp2$betafrailgap-(-1))^2, na.rm=TRUE)
msqer31 <- mean((exp3$betafrailgap-(-1))^2, na.rm=TRUE)
msqer41 <- mean((exp4$betafrailgap-(-1))^2, na.rm=TRUE)
msqer51 <- mean((exp5$betafrailgap-(-1))^2, na.rm=TRUE)
msqer61 <- mean((exp6$betafrailgap-(-1))^2, na.rm=TRUE)
msqer71 <- mean((exp7$betafrailgap-(-1))^2, na.rm=TRUE)
msqer81 <- mean((exp8$betafrailgap-(-1))^2, na.rm=TRUE)
msqer91 <- mean((exp9$betafrailgap-(-1))^2, na.rm=TRUE)


msqer12 <- mean((exp1$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer22 <- mean((exp2$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer32 <- mean((exp3$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer42 <- mean((exp4$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer52 <- mean((exp5$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer62 <- mean((exp6$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer72 <- mean((exp7$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer82 <- mean((exp8$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer92 <- mean((exp9$betafrailelapse-(-1))^2, na.rm=TRUE)


msqer13 <- mean((exp1$betaag-(-1))^2, na.rm=TRUE)
msqer23 <- mean((exp2$betaag-(-1))^2, na.rm=TRUE)
msqer33 <- mean((exp3$betaag-(-1))^2, na.rm=TRUE)
msqer43 <- mean((exp4$betaag-(-1))^2, na.rm=TRUE)
msqer53 <- mean((exp5$betaag-(-1))^2, na.rm=TRUE)
msqer63 <- mean((exp6$betaag-(-1))^2, na.rm=TRUE)
msqer73 <- mean((exp7$betaag-(-1))^2, na.rm=TRUE)
msqer83 <- mean((exp8$betaag-(-1))^2, na.rm=TRUE)
msqer93 <- mean((exp9$betaag-(-1))^2, na.rm=TRUE)


msqer14 <- mean((exp1$betacondgap-(-1))^2, na.rm=TRUE)
msqer24 <- mean((exp2$betacondgap-(-1))^2, na.rm=TRUE)
msqer34 <- mean((exp3$betacondgap-(-1))^2, na.rm=TRUE)
msqer44 <- mean((exp4$betacondgap-(-1))^2, na.rm=TRUE)
msqer54 <- mean((exp5$betacondgap-(-1))^2, na.rm=TRUE)
msqer64 <- mean((exp6$betacondgap-(-1))^2, na.rm=TRUE)
msqer74 <- mean((exp7$betacondgap-(-1))^2, na.rm=TRUE)
msqer84 <- mean((exp8$betacondgap-(-1))^2, na.rm=TRUE)
msqer94 <- mean((exp9$betacondgap-(-1))^2, na.rm=TRUE)


msqer15 <- mean((exp1$betacondelapse-(-1))^2, na.rm=TRUE)
msqer25 <- mean((exp2$betacondelapse-(-1))^2, na.rm=TRUE)
msqer35 <- mean((exp3$betacondelapse-(-1))^2, na.rm=TRUE)
msqer45 <- mean((exp4$betacondelapse-(-1))^2, na.rm=TRUE)
msqer55 <- mean((exp5$betacondelapse-(-1))^2, na.rm=TRUE)
msqer65 <- mean((exp6$betacondelapse-(-1))^2, na.rm=TRUE)
msqer75 <- mean((exp7$betacondelapse-(-1))^2, na.rm=TRUE)
msqer85 <- mean((exp8$betacondelapse-(-1))^2, na.rm=TRUE)
msqer95 <- mean((exp9$betacondelapse-(-1))^2, na.rm=TRUE)


msqcondfrail <- c(msqer11, msqer21, msqer31, msqer41,msqer51,msqer61,msqer71, msqer81, msqer91)
msqag <-  c(msqer13, msqer23, msqer33, msqer43,msqer53,msqer63,msqer73, msqer83, msqer93)
msqfrailelapse <- c(msqer12, msqer22, msqer32, msqer42,msqer52,msqer62,msqer72, msqer82, msqer92)
msqcondgap <- c(msqer14, msqer24, msqer34, msqer44,msqer54,msqer64,msqer74, msqer84, msqer94)
msqcondelapse <- c(msqer15, msqer25, msqer35, msqer45,msqer55,msqer65,msqer75, msqer85, msqer95)

sample <- c(2,3,4,5,6,7,8,9,10)

yvalues <- c(msqcondfrail,msqcondgap,msqcondelapse,msqag,msqfrailelapse)
xvalues <- c(rep(sample,5))
names <- c(rep("Conditional Frailty",9), rep("Conditional, Gap",9), rep("Conditional, Elapse",9), rep("Andersen-Gill",9), rep("Frailty, Elapse",9))
ednohet <- data.frame(msqerror=yvalues,sample=xvalues,Model=names)

exp1 <- read.table(file="exp_noed_het_maxev2.txt", header=TRUE)
exp2 <- read.table(file="exp_noed_het_maxev3.txt", header=TRUE)
exp3 <- read.table(file="exp_noed_het_maxev4.txt", header=TRUE)
exp4 <- read.table(file="exp_noed_het_maxev5.txt", header=TRUE)
exp5 <- read.table(file="exp_noed_het_maxev6.txt", header=TRUE)
exp6 <- read.table(file="exp_noed_het_maxev7.txt", header=TRUE)
exp7 <- read.table(file="exp_noed_het_maxev8.txt", header=TRUE)
exp8 <- read.table(file="exp_noed_het_maxev9.txt", header=TRUE)
exp9 <- read.table(file="exp_noed_het_maxev10.txt", header=TRUE)

msqer11 <- mean((exp1$betafrailgap-(-1))^2, na.rm=TRUE)
msqer21 <- mean((exp2$betafrailgap-(-1))^2, na.rm=TRUE)
msqer31 <- mean((exp3$betafrailgap-(-1))^2, na.rm=TRUE)
msqer41 <- mean((exp4$betafrailgap-(-1))^2, na.rm=TRUE)
msqer51 <- mean((exp5$betafrailgap-(-1))^2, na.rm=TRUE)
msqer61 <- mean((exp6$betafrailgap-(-1))^2, na.rm=TRUE)
msqer71 <- mean((exp7$betafrailgap-(-1))^2, na.rm=TRUE)
msqer81 <- mean((exp8$betafrailgap-(-1))^2, na.rm=TRUE)
msqer91 <- mean((exp9$betafrailgap-(-1))^2, na.rm=TRUE)


msqer12 <- mean((exp1$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer22 <- mean((exp2$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer32 <- mean((exp3$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer42 <- mean((exp4$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer52 <- mean((exp5$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer62 <- mean((exp6$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer72 <- mean((exp7$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer82 <- mean((exp8$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer92 <- mean((exp9$betafrailelapse-(-1))^2, na.rm=TRUE)


msqer13 <- mean((exp1$betaag-(-1))^2, na.rm=TRUE)
msqer23 <- mean((exp2$betaag-(-1))^2, na.rm=TRUE)
msqer33 <- mean((exp3$betaag-(-1))^2, na.rm=TRUE)
msqer43 <- mean((exp4$betaag-(-1))^2, na.rm=TRUE)
msqer53 <- mean((exp5$betaag-(-1))^2, na.rm=TRUE)
msqer63 <- mean((exp6$betaag-(-1))^2, na.rm=TRUE)
msqer73 <- mean((exp7$betaag-(-1))^2, na.rm=TRUE)
msqer83 <- mean((exp8$betaag-(-1))^2, na.rm=TRUE)
msqer93 <- mean((exp9$betaag-(-1))^2, na.rm=TRUE)


msqer14 <- mean((exp1$betacondgap-(-1))^2, na.rm=TRUE)
msqer24 <- mean((exp2$betacondgap-(-1))^2, na.rm=TRUE)
msqer34 <- mean((exp3$betacondgap-(-1))^2, na.rm=TRUE)
msqer44 <- mean((exp4$betacondgap-(-1))^2, na.rm=TRUE)
msqer54 <- mean((exp5$betacondgap-(-1))^2, na.rm=TRUE)
msqer64 <- mean((exp6$betacondgap-(-1))^2, na.rm=TRUE)
msqer74 <- mean((exp7$betacondgap-(-1))^2, na.rm=TRUE)
msqer84 <- mean((exp8$betacondgap-(-1))^2, na.rm=TRUE)
msqer94 <- mean((exp9$betacondgap-(-1))^2, na.rm=TRUE)


msqer15 <- mean((exp1$betacondelapse-(-1))^2, na.rm=TRUE)
msqer25 <- mean((exp2$betacondelapse-(-1))^2, na.rm=TRUE)
msqer35 <- mean((exp3$betacondelapse-(-1))^2, na.rm=TRUE)
msqer45 <- mean((exp4$betacondelapse-(-1))^2, na.rm=TRUE)
msqer55 <- mean((exp5$betacondelapse-(-1))^2, na.rm=TRUE)
msqer65 <- mean((exp6$betacondelapse-(-1))^2, na.rm=TRUE)
msqer75 <- mean((exp7$betacondelapse-(-1))^2, na.rm=TRUE)
msqer85 <- mean((exp8$betacondelapse-(-1))^2, na.rm=TRUE)
msqer95 <- mean((exp9$betacondelapse-(-1))^2, na.rm=TRUE)


msqcondfrail <- c(msqer11, msqer21, msqer31, msqer41,msqer51,msqer61,msqer71, msqer81, msqer91)
msqag <-  c(msqer13, msqer23, msqer33, msqer43,msqer53,msqer63,msqer73, msqer83, msqer93)
msqfrailelapse <- c(msqer12, msqer22, msqer32, msqer42,msqer52,msqer62,msqer72, msqer82, msqer92)
msqcondgap <- c(msqer14, msqer24, msqer34, msqer44,msqer54,msqer64,msqer74, msqer84, msqer94)
msqcondelapse <- c(msqer15, msqer25, msqer35, msqer45,msqer55,msqer65,msqer75, msqer85, msqer95)
      
sample <- c(2,3,4,5,6,7,8,9,10)

yvalues <- c(msqcondfrail,msqcondgap,msqcondelapse,msqag,msqfrailelapse)
xvalues <- c(rep(sample,5))
names <- c(rep("Conditional Frailty",9), rep("Conditional, Gap",9), rep("Conditional, Elapse",9), rep("Andersen-Gill",9), rep("Frailty, Elapse",9))
noedhet <- data.frame(msqerror=yvalues,sample=xvalues,Model=names)

exp1 <- read.table(file="exp_ed_het_maxev2.txt", header=TRUE)
exp2 <- read.table(file="exp_ed_het_maxev3.txt", header=TRUE)
exp3 <- read.table(file="exp_ed_het_maxev4.txt", header=TRUE)
exp4 <- read.table(file="exp_ed_het_maxev5.txt", header=TRUE)
exp5 <- read.table(file="exp_ed_het_maxev6.txt", header=TRUE)
exp6 <- read.table(file="exp_ed_het_maxev7.txt", header=TRUE)
exp7 <- read.table(file="exp_ed_het_maxev8.txt", header=TRUE)
exp8 <- read.table(file="exp_ed_het_maxev9.txt", header=TRUE)
exp9 <- read.table(file="exp_ed_het_maxev10.txt", header=TRUE)

msqer11 <- mean((exp1$betafrailgap-(-1))^2, na.rm=TRUE)
msqer21 <- mean((exp2$betafrailgap-(-1))^2, na.rm=TRUE)
msqer31 <- mean((exp3$betafrailgap-(-1))^2, na.rm=TRUE)
msqer41 <- mean((exp4$betafrailgap-(-1))^2, na.rm=TRUE)
msqer51 <- mean((exp5$betafrailgap-(-1))^2, na.rm=TRUE)
msqer61 <- mean((exp6$betafrailgap-(-1))^2, na.rm=TRUE)
msqer71 <- mean((exp7$betafrailgap-(-1))^2, na.rm=TRUE)
msqer81 <- mean((exp8$betafrailgap-(-1))^2, na.rm=TRUE)
msqer91 <- mean((exp9$betafrailgap-(-1))^2, na.rm=TRUE)


msqer12 <- mean((exp1$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer22 <- mean((exp2$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer32 <- mean((exp3$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer42 <- mean((exp4$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer52 <- mean((exp5$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer62 <- mean((exp6$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer72 <- mean((exp7$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer82 <- mean((exp8$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer92 <- mean((exp9$betafrailelapse-(-1))^2, na.rm=TRUE)


msqer13 <- mean((exp1$betaag-(-1))^2, na.rm=TRUE)
msqer23 <- mean((exp2$betaag-(-1))^2, na.rm=TRUE)
msqer33 <- mean((exp3$betaag-(-1))^2, na.rm=TRUE)
msqer43 <- mean((exp4$betaag-(-1))^2, na.rm=TRUE)
msqer53 <- mean((exp5$betaag-(-1))^2, na.rm=TRUE)
msqer63 <- mean((exp6$betaag-(-1))^2, na.rm=TRUE)
msqer73 <- mean((exp7$betaag-(-1))^2, na.rm=TRUE)
msqer83 <- mean((exp8$betaag-(-1))^2, na.rm=TRUE)
msqer93 <- mean((exp9$betaag-(-1))^2, na.rm=TRUE)


msqer14 <- mean((exp1$betacondgap-(-1))^2, na.rm=TRUE)
msqer24 <- mean((exp2$betacondgap-(-1))^2, na.rm=TRUE)
msqer34 <- mean((exp3$betacondgap-(-1))^2, na.rm=TRUE)
msqer44 <- mean((exp4$betacondgap-(-1))^2, na.rm=TRUE)
msqer54 <- mean((exp5$betacondgap-(-1))^2, na.rm=TRUE)
msqer64 <- mean((exp6$betacondgap-(-1))^2, na.rm=TRUE)
msqer74 <- mean((exp7$betacondgap-(-1))^2, na.rm=TRUE)
msqer84 <- mean((exp8$betacondgap-(-1))^2, na.rm=TRUE)
msqer94 <- mean((exp9$betacondgap-(-1))^2, na.rm=TRUE)


msqer15 <- mean((exp1$betacondelapse-(-1))^2, na.rm=TRUE)
msqer25 <- mean((exp2$betacondelapse-(-1))^2, na.rm=TRUE)
msqer35 <- mean((exp3$betacondelapse-(-1))^2, na.rm=TRUE)
msqer45 <- mean((exp4$betacondelapse-(-1))^2, na.rm=TRUE)
msqer55 <- mean((exp5$betacondelapse-(-1))^2, na.rm=TRUE)
msqer65 <- mean((exp6$betacondelapse-(-1))^2, na.rm=TRUE)
msqer75 <- mean((exp7$betacondelapse-(-1))^2, na.rm=TRUE)
msqer85 <- mean((exp8$betacondelapse-(-1))^2, na.rm=TRUE)
msqer95 <- mean((exp9$betacondelapse-(-1))^2, na.rm=TRUE)


msqcondfrail <- c(msqer11, msqer21, msqer31, msqer41,msqer51,msqer61,msqer71, msqer81, msqer91)
msqag <-  c(msqer13, msqer23, msqer33, msqer43,msqer53,msqer63,msqer73, msqer83, msqer93)
msqfrailelapse <- c(msqer12, msqer22, msqer32, msqer42,msqer52,msqer62,msqer72, msqer82, msqer92)
msqcondgap <- c(msqer14, msqer24, msqer34, msqer44,msqer54,msqer64,msqer74, msqer84, msqer94)
msqcondelapse <- c(msqer15, msqer25, msqer35, msqer45,msqer55,msqer65,msqer75, msqer85, msqer95)
      
sample <- c(2,3,4,5,6,7,8,9,10)


yvalues <- c(msqcondfrail,msqcondgap,msqcondelapse,msqag,msqfrailelapse)
xvalues <- c(rep(sample,5))
names <- c(rep("Conditional Frailty",9), rep("Conditional, Gap",9), rep("Conditional, Elapse",9), rep("Andersen-Gill",9), rep("Frailty, Elapse",9))
edhet <- data.frame(msqerror=yvalues,sample=xvalues,Model=names)


msqplot1 <- qplot(sample,msqerror,data=noednohet,geom="line",colour=Model,linetype=Model) + scale_colour_brewer(pal="Dark2") + labs(y="Mean Squared Error", x="Maximum Events (Cases=100)") + opts(legend.position="none") + scale_linetype_manual("Model",values=c(1:2,4:6)) + opts(title="No Event Dependence, No Heterogeneity")

msqplot2 <- qplot(sample,msqerror,data=ednohet,geom="line",colour=Model,linetype=Model) + scale_colour_brewer(pal="Dark2") + labs(y="Mean Squared Error", x="Maximum Events (Cases=100)", main="Event Dependence, No Heterogeneity") + ylim(-.01,.6) + opts(legend.position=c(.7,.7)) + scale_linetype_manual("Model",values=c(1:2,4:6))  + opts(title="Event Dependence, No Heterogeneity")

msqplot3 <- qplot(sample,msqerror,data=noedhet,geom="line",colour=Model,linetype=Model) + scale_colour_brewer(pal="Dark2") + labs(y="Mean Squared Error", x="Maximum Events (Cases=100)", main="No Event Dependence, Heterogeneity") + opts(legend.position="none")  + scale_linetype_manual("Model",values=c(1:2,4:6))  + opts(title="No Event Dependence, Heterogeneity")

msqplot4 <- qplot(sample,msqerror,data=edhet,geom="line",colour=Model,linetype=Model) + scale_colour_brewer(pal="Dark2") + labs(y="Mean Squared Error", x="Maximum Events (Cases=100)", main="Event Dependence, Heterogeneity") + opts(legend.position="none")  + scale_linetype_manual("Model",values=c(1:2,4:6))  + opts(title="Event Dependence, Heterogeneity")


## B & W


msqplot1 <- qplot(sample,msqerror,data=noednohet,geom="line",linetype=Model) + theme_bw()+  labs(y="Mean Squared Error", x="Maximum Events (Cases=100)") + opts(legend.position=c(.7,.7)) + scale_linetype_manual("Model",values=c(1,2,4,3,5)) + opts(title="No Event Dependence, No Heterogeneity") +ylim(-.01,.62)

msqplot2 <- qplot(sample,msqerror,data=ednohet,geom="line",linetype=Model) + theme_bw() +  labs(y="Mean Squared Error", x="Maximum Events (Cases=100)")  + opts(legend.position="none") + scale_linetype_manual("Model",values=c(1,2,4,3,5))  + opts(title="Event Dependence, No Heterogeneity")  +ylim(-.01,.62)

msqplot3 <- qplot(sample,msqerror,data=noedhet,geom="line",linetype=Model) + theme_bw() +  labs(y="Mean Squared Error", x="Maximum Events (Cases=100)") + opts(legend.position="none")  + scale_linetype_manual("Model",values=c(1,2,4,3,5))  + opts(title="No Event Dependence, Heterogeneity")   +ylim(-.01,.62)

msqplot4 <- qplot(sample,msqerror,data=edhet,geom="line",linetype=Model) + theme_bw()+  labs(y="Mean Squared Error", x="Maximum Events (Cases=100)") + opts(legend.position="none")  + scale_linetype_manual("Model",values=c(1,2,4,3,5))  + opts(title="Event Dependence, Heterogeneity")  +ylim(-.01,.62)



pdf("mse_maxev_panel.pdf",height=8.5,width=11)
grid.newpage()
pushViewport(viewport(layout = grid.layout(2,2)))
vplayout <- function(x,y)
viewport(layout.pos.row=x,layout.pos.col=y)
print(msqplot1,vp=vplayout(1,1))
print(msqplot2,vp=vplayout(1,2))
print(msqplot3,vp=vplayout(2,1))
print(msqplot4,vp=vplayout(2,2))
dev.off()


exp1 <- read.table(file="exp_noed_nohet_samplesize_ver210.txt", header=TRUE)
exp2 <- read.table(file="exp_noed_nohet_samplesize_ver230.txt", header=TRUE)
exp3 <- read.table(file="exp_noed_nohet_samplesize_ver250.txt", header=TRUE)
exp4 <- read.table(file="exp_noed_nohet_samplesize_ver270.txt", header=TRUE)
exp5 <- read.table(file="exp_noed_nohet_samplesize_ver290.txt", header=TRUE)
exp6 <- read.table(file="exp_noed_nohet_samplesize_ver2110.txt", header=TRUE)
exp7 <- read.table(file="exp_noed_nohet_samplesize_ver2130.txt", header=TRUE)
exp8 <- read.table(file="exp_noed_nohet_samplesize_ver2150.txt", header=TRUE)
exp9 <- read.table(file="exp_noed_nohet_samplesize_ver2170.txt", header=TRUE)
exp10 <- read.table(file="exp_noed_nohet_samplesize_ver2190.txt", header=TRUE)
exp11 <- read.table(file="exp_noed_nohet_samplesize_ver2210.txt", header=TRUE)
exp12 <- read.table(file="exp_noed_nohet_samplesize_ver2230.txt", header=TRUE)
exp13 <- read.table(file="exp_noed_nohet_samplesize_ver2250.txt", header=TRUE)
exp14 <- read.table(file="exp_noed_nohet_samplesize_ver2270.txt", header=TRUE)
exp15 <- read.table(file="exp_noed_nohet_samplesize_ver2290.txt", header=TRUE)
exp16 <- read.table(file="exp_noed_nohet_samplesize_ver2310.txt", header=TRUE)
exp17 <- read.table(file="exp_noed_nohet_samplesize_ver2330.txt", header=TRUE)


msqer11 <- mean((exp1$betafrailgap-(-1))^2, na.rm=TRUE)
msqer21 <- mean((exp2$betafrailgap-(-1))^2, na.rm=TRUE)
msqer31 <- mean((exp3$betafrailgap-(-1))^2, na.rm=TRUE)
msqer41 <- mean((exp4$betafrailgap-(-1))^2, na.rm=TRUE)
msqer51 <- mean((exp5$betafrailgap-(-1))^2, na.rm=TRUE)
msqer61 <- mean((exp6$betafrailgap-(-1))^2, na.rm=TRUE)
msqer71 <- mean((exp7$betafrailgap-(-1))^2, na.rm=TRUE)
msqer81 <- mean((exp8$betafrailgap-(-1))^2, na.rm=TRUE)
msqer91 <- mean((exp9$betafrailgap-(-1))^2, na.rm=TRUE)
msqer101 <- mean((exp10$betafrailgap-(-1))^2, na.rm=TRUE)
msqer111 <- mean((exp11$betafrailgap-(-1))^2, na.rm=TRUE)
msqer121 <- mean((exp12$betafrailgap-(-1))^2, na.rm=TRUE)
msqer131 <- mean((exp13$betafrailgap-(-1))^2, na.rm=TRUE)
msqer141 <- mean((exp14$betafrailgap-(-1))^2, na.rm=TRUE)
msqer151 <- mean((exp15$betafrailgap-(-1))^2, na.rm=TRUE)
msqer161 <- mean((exp16$betafrailgap-(-1))^2, na.rm=TRUE)
msqer171 <- mean((exp17$betafrailgap-(-1))^2, na.rm=TRUE)

msqer12 <- mean((exp1$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer22 <- mean((exp2$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer32 <- mean((exp3$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer42 <- mean((exp4$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer52 <- mean((exp5$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer62 <- mean((exp6$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer72 <- mean((exp7$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer82 <- mean((exp8$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer92 <- mean((exp9$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer102 <- mean((exp10$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer112 <- mean((exp11$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer122 <- mean((exp12$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer132 <- mean((exp13$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer142 <- mean((exp14$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer152 <- mean((exp15$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer162 <- mean((exp16$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer172 <- mean((exp17$betafrailelapse-(-1))^2, na.rm=TRUE)

msqer13 <- mean((exp1$betaag-(-1))^2, na.rm=TRUE)
msqer23 <- mean((exp2$betaag-(-1))^2, na.rm=TRUE)
msqer33 <- mean((exp3$betaag-(-1))^2, na.rm=TRUE)
msqer43 <- mean((exp4$betaag-(-1))^2, na.rm=TRUE)
msqer53 <- mean((exp5$betaag-(-1))^2, na.rm=TRUE)
msqer63 <- mean((exp6$betaag-(-1))^2, na.rm=TRUE)
msqer73 <- mean((exp7$betaag-(-1))^2, na.rm=TRUE)
msqer83 <- mean((exp8$betaag-(-1))^2, na.rm=TRUE)
msqer93 <- mean((exp9$betaag-(-1))^2, na.rm=TRUE)
msqer103 <- mean((exp10$betaag-(-1))^2, na.rm=TRUE)
msqer113 <- mean((exp11$betaag-(-1))^2, na.rm=TRUE)
msqer123 <- mean((exp12$betaag-(-1))^2, na.rm=TRUE)
msqer133 <- mean((exp13$betaag-(-1))^2, na.rm=TRUE)
msqer143 <- mean((exp14$betaag-(-1))^2, na.rm=TRUE)
msqer153 <- mean((exp15$betaag-(-1))^2, na.rm=TRUE)
msqer163 <- mean((exp16$betaag-(-1))^2, na.rm=TRUE)
msqer173 <- mean((exp17$betaag-(-1))^2, na.rm=TRUE)

msqer14 <- mean((exp1$betacondgap-(-1))^2, na.rm=TRUE)
msqer24 <- mean((exp2$betacondgap-(-1))^2, na.rm=TRUE)
msqer34 <- mean((exp3$betacondgap-(-1))^2, na.rm=TRUE)
msqer44 <- mean((exp4$betacondgap-(-1))^2, na.rm=TRUE)
msqer54 <- mean((exp5$betacondgap-(-1))^2, na.rm=TRUE)
msqer64 <- mean((exp6$betacondgap-(-1))^2, na.rm=TRUE)
msqer74 <- mean((exp7$betacondgap-(-1))^2, na.rm=TRUE)
msqer84 <- mean((exp8$betacondgap-(-1))^2, na.rm=TRUE)
msqer94 <- mean((exp9$betacondgap-(-1))^2, na.rm=TRUE)
msqer104 <- mean((exp10$betacondgap-(-1))^2, na.rm=TRUE)
msqer114 <- mean((exp11$betacondgap-(-1))^2, na.rm=TRUE)
msqer124 <- mean((exp12$betacondgap-(-1))^2, na.rm=TRUE)
msqer134 <- mean((exp13$betacondgap-(-1))^2, na.rm=TRUE)
msqer144 <- mean((exp14$betacondgap-(-1))^2, na.rm=TRUE)
msqer154 <- mean((exp15$betacondgap-(-1))^2, na.rm=TRUE)
msqer164 <- mean((exp16$betacondgap-(-1))^2, na.rm=TRUE)
msqer174 <- mean((exp17$betacondgap-(-1))^2, na.rm=TRUE)

msqer15 <- mean((exp1$betacondelapse-(-1))^2, na.rm=TRUE)
msqer25 <- mean((exp2$betacondelapse-(-1))^2, na.rm=TRUE)
msqer35 <- mean((exp3$betacondelapse-(-1))^2, na.rm=TRUE)
msqer45 <- mean((exp4$betacondelapse-(-1))^2, na.rm=TRUE)
msqer55 <- mean((exp5$betacondelapse-(-1))^2, na.rm=TRUE)
msqer65 <- mean((exp6$betacondelapse-(-1))^2, na.rm=TRUE)
msqer75 <- mean((exp7$betacondelapse-(-1))^2, na.rm=TRUE)
msqer85 <- mean((exp8$betacondelapse-(-1))^2, na.rm=TRUE)
msqer95 <- mean((exp9$betacondelapse-(-1))^2, na.rm=TRUE)
msqer105 <- mean((exp10$betacondelapse-(-1))^2, na.rm=TRUE)
msqer115 <- mean((exp11$betacondelapse-(-1))^2, na.rm=TRUE)
msqer125 <- mean((exp12$betacondelapse-(-1))^2, na.rm=TRUE)
msqer135 <- mean((exp13$betacondelapse-(-1))^2, na.rm=TRUE)
msqer145 <- mean((exp14$betacondelapse-(-1))^2, na.rm=TRUE)
msqer155 <- mean((exp15$betacondelapse-(-1))^2, na.rm=TRUE)
msqer165 <- mean((exp16$betacondelapse-(-1))^2, na.rm=TRUE)
msqer175 <- mean((exp17$betacondelapse-(-1))^2, na.rm=TRUE)

msqcondfrail <- c(msqer11, msqer21, msqer31, msqer41,msqer51,msqer61,msqer71, msqer81, msqer91, msqer101,msqer111, msqer121, msqer131, msqer141,msqer151,msqer161,msqer171)
msqag <-  c(msqer13, msqer23, msqer33, msqer43,msqer53,msqer63,msqer73, msqer83, msqer93,msqer103,msqer113, msqer123, msqer133, msqer143,msqer153,msqer163,msqer173)
msqfrailelapse <- c(msqer12, msqer22, msqer32, msqer42,msqer52,msqer62,msqer72, msqer82, msqer92,msqer102, msqer112, msqer122, msqer132, msqer142,msqer152,msqer162,msqer172)
msqcondgap <- c(msqer14, msqer24, msqer34, msqer44,msqer54,msqer64,msqer74, msqer84, msqer94, msqer104, msqer114, msqer124, msqer134, msqer144,msqer154,msqer164,msqer174)
msqcondelapse <- c(msqer15, msqer25, msqer35, msqer45,msqer55,msqer65,msqer75, msqer85, msqer95, msqer105, msqer115, msqer125, msqer135, msqer145,msqer155,msqer165,msqer175)
      
sample <- c(10,30,50,70,90,110,130,150,170,190,210,230,250,270,290,310,330)

yvalues <- c(msqcondfrail,msqcondgap,msqcondelapse,msqag,msqfrailelapse)
xvalues <- c(rep(sample,5))
names <- c(rep("Conditional Frailty",17), rep("Conditional, Gap",17), rep("Conditional, Elapse",17), rep("Andersen-Gill",17), rep("Frailty, Elapse",17))
noednohet <- data.frame(msqerror=yvalues,sample=xvalues,Model=names)


exp1 <- read.table(file="exp_ed_nohet_samplesize_ver210.txt", header=TRUE)
exp2 <- read.table(file="exp_ed_nohet_samplesize_ver230.txt", header=TRUE)
exp3 <- read.table(file="exp_ed_nohet_samplesize_ver250.txt", header=TRUE)
exp4 <- read.table(file="exp_ed_nohet_samplesize_ver270.txt", header=TRUE)
exp5 <- read.table(file="exp_ed_nohet_samplesize_ver290.txt", header=TRUE)
exp6 <- read.table(file="exp_ed_nohet_samplesize_ver2110.txt", header=TRUE)
exp7 <- read.table(file="exp_ed_nohet_samplesize_ver2130.txt", header=TRUE)
exp8 <- read.table(file="exp_ed_nohet_samplesize_ver2150.txt", header=TRUE)
exp9 <- read.table(file="exp_ed_nohet_samplesize_ver2170.txt", header=TRUE)
exp10 <- read.table(file="exp_ed_nohet_samplesize_ver2190.txt", header=TRUE)
exp11 <- read.table(file="exp_ed_nohet_samplesize_ver2210.txt", header=TRUE)
exp12 <- read.table(file="exp_ed_nohet_samplesize_ver2230.txt", header=TRUE)
exp13 <- read.table(file="exp_ed_nohet_samplesize_ver2250.txt", header=TRUE)
exp14 <- read.table(file="exp_ed_nohet_samplesize_ver2270.txt", header=TRUE)
exp15 <- read.table(file="exp_ed_nohet_samplesize_ver2290.txt", header=TRUE)
exp16 <- read.table(file="exp_ed_nohet_samplesize_ver2310.txt", header=TRUE)
exp17 <- read.table(file="exp_ed_nohet_samplesize_ver2330.txt", header=TRUE)


msqer11 <- mean((exp1$betafrailgap-(-1))^2, na.rm=TRUE)
msqer21 <- mean((exp2$betafrailgap-(-1))^2, na.rm=TRUE)
msqer31 <- mean((exp3$betafrailgap-(-1))^2, na.rm=TRUE)
msqer41 <- mean((exp4$betafrailgap-(-1))^2, na.rm=TRUE)
msqer51 <- mean((exp5$betafrailgap-(-1))^2, na.rm=TRUE)
msqer61 <- mean((exp6$betafrailgap-(-1))^2, na.rm=TRUE)
msqer71 <- mean((exp7$betafrailgap-(-1))^2, na.rm=TRUE)
msqer81 <- mean((exp8$betafrailgap-(-1))^2, na.rm=TRUE)
msqer91 <- mean((exp9$betafrailgap-(-1))^2, na.rm=TRUE)
msqer101 <- mean((exp10$betafrailgap-(-1))^2, na.rm=TRUE)
msqer111 <- mean((exp11$betafrailgap-(-1))^2, na.rm=TRUE)
msqer121 <- mean((exp12$betafrailgap-(-1))^2, na.rm=TRUE)
msqer131 <- mean((exp13$betafrailgap-(-1))^2, na.rm=TRUE)
msqer141 <- mean((exp14$betafrailgap-(-1))^2, na.rm=TRUE)
msqer151 <- mean((exp15$betafrailgap-(-1))^2, na.rm=TRUE)
msqer161 <- mean((exp16$betafrailgap-(-1))^2, na.rm=TRUE)
msqer171 <- mean((exp17$betafrailgap-(-1))^2, na.rm=TRUE)

msqer12 <- mean((exp1$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer22 <- mean((exp2$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer32 <- mean((exp3$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer42 <- mean((exp4$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer52 <- mean((exp5$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer62 <- mean((exp6$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer72 <- mean((exp7$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer82 <- mean((exp8$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer92 <- mean((exp9$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer102 <- mean((exp10$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer112 <- mean((exp11$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer122 <- mean((exp12$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer132 <- mean((exp13$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer142 <- mean((exp14$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer152 <- mean((exp15$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer162 <- mean((exp16$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer172 <- mean((exp17$betafrailelapse-(-1))^2, na.rm=TRUE)

msqer13 <- mean((exp1$betaag-(-1))^2, na.rm=TRUE)
msqer23 <- mean((exp2$betaag-(-1))^2, na.rm=TRUE)
msqer33 <- mean((exp3$betaag-(-1))^2, na.rm=TRUE)
msqer43 <- mean((exp4$betaag-(-1))^2, na.rm=TRUE)
msqer53 <- mean((exp5$betaag-(-1))^2, na.rm=TRUE)
msqer63 <- mean((exp6$betaag-(-1))^2, na.rm=TRUE)
msqer73 <- mean((exp7$betaag-(-1))^2, na.rm=TRUE)
msqer83 <- mean((exp8$betaag-(-1))^2, na.rm=TRUE)
msqer93 <- mean((exp9$betaag-(-1))^2, na.rm=TRUE)
msqer103 <- mean((exp10$betaag-(-1))^2, na.rm=TRUE)
msqer113 <- mean((exp11$betaag-(-1))^2, na.rm=TRUE)
msqer123 <- mean((exp12$betaag-(-1))^2, na.rm=TRUE)
msqer133 <- mean((exp13$betaag-(-1))^2, na.rm=TRUE)
msqer143 <- mean((exp14$betaag-(-1))^2, na.rm=TRUE)
msqer153 <- mean((exp15$betaag-(-1))^2, na.rm=TRUE)
msqer163 <- mean((exp16$betaag-(-1))^2, na.rm=TRUE)
msqer173 <- mean((exp17$betaag-(-1))^2, na.rm=TRUE)

msqer14 <- mean((exp1$betacondgap-(-1))^2, na.rm=TRUE)
msqer24 <- mean((exp2$betacondgap-(-1))^2, na.rm=TRUE)
msqer34 <- mean((exp3$betacondgap-(-1))^2, na.rm=TRUE)
msqer44 <- mean((exp4$betacondgap-(-1))^2, na.rm=TRUE)
msqer54 <- mean((exp5$betacondgap-(-1))^2, na.rm=TRUE)
msqer64 <- mean((exp6$betacondgap-(-1))^2, na.rm=TRUE)
msqer74 <- mean((exp7$betacondgap-(-1))^2, na.rm=TRUE)
msqer84 <- mean((exp8$betacondgap-(-1))^2, na.rm=TRUE)
msqer94 <- mean((exp9$betacondgap-(-1))^2, na.rm=TRUE)
msqer104 <- mean((exp10$betacondgap-(-1))^2, na.rm=TRUE)
msqer114 <- mean((exp11$betacondgap-(-1))^2, na.rm=TRUE)
msqer124 <- mean((exp12$betacondgap-(-1))^2, na.rm=TRUE)
msqer134 <- mean((exp13$betacondgap-(-1))^2, na.rm=TRUE)
msqer144 <- mean((exp14$betacondgap-(-1))^2, na.rm=TRUE)
msqer154 <- mean((exp15$betacondgap-(-1))^2, na.rm=TRUE)
msqer164 <- mean((exp16$betacondgap-(-1))^2, na.rm=TRUE)
msqer174 <- mean((exp17$betacondgap-(-1))^2, na.rm=TRUE)

msqer15 <- mean((exp1$betacondelapse-(-1))^2, na.rm=TRUE)
msqer25 <- mean((exp2$betacondelapse-(-1))^2, na.rm=TRUE)
msqer35 <- mean((exp3$betacondelapse-(-1))^2, na.rm=TRUE)
msqer45 <- mean((exp4$betacondelapse-(-1))^2, na.rm=TRUE)
msqer55 <- mean((exp5$betacondelapse-(-1))^2, na.rm=TRUE)
msqer65 <- mean((exp6$betacondelapse-(-1))^2, na.rm=TRUE)
msqer75 <- mean((exp7$betacondelapse-(-1))^2, na.rm=TRUE)
msqer85 <- mean((exp8$betacondelapse-(-1))^2, na.rm=TRUE)
msqer95 <- mean((exp9$betacondelapse-(-1))^2, na.rm=TRUE)
msqer105 <- mean((exp10$betacondelapse-(-1))^2, na.rm=TRUE)
msqer115 <- mean((exp11$betacondelapse-(-1))^2, na.rm=TRUE)
msqer125 <- mean((exp12$betacondelapse-(-1))^2, na.rm=TRUE)
msqer135 <- mean((exp13$betacondelapse-(-1))^2, na.rm=TRUE)
msqer145 <- mean((exp14$betacondelapse-(-1))^2, na.rm=TRUE)
msqer155 <- mean((exp15$betacondelapse-(-1))^2, na.rm=TRUE)
msqer165 <- mean((exp16$betacondelapse-(-1))^2, na.rm=TRUE)
msqer175 <- mean((exp17$betacondelapse-(-1))^2, na.rm=TRUE)

msqcondfrail <- c(msqer11, msqer21, msqer31, msqer41,msqer51,msqer61,msqer71, msqer81, msqer91, msqer101,msqer111, msqer121, msqer131, msqer141,msqer151,msqer161,msqer171)
msqag <-  c(msqer13, msqer23, msqer33, msqer43,msqer53,msqer63,msqer73, msqer83, msqer93,msqer103,msqer113, msqer123, msqer133, msqer143,msqer153,msqer163,msqer173)
msqfrailelapse <- c(msqer12, msqer22, msqer32, msqer42,msqer52,msqer62,msqer72, msqer82, msqer92,msqer102, msqer112, msqer122, msqer132, msqer142,msqer152,msqer162,msqer172)
msqcondgap <- c(msqer14, msqer24, msqer34, msqer44,msqer54,msqer64,msqer74, msqer84, msqer94, msqer104, msqer114, msqer124, msqer134, msqer144,msqer154,msqer164,msqer174)
msqcondelapse <- c(msqer15, msqer25, msqer35, msqer45,msqer55,msqer65,msqer75, msqer85, msqer95, msqer105, msqer115, msqer125, msqer135, msqer145,msqer155,msqer165,msqer175)
      
sample <- c(10,30,50,70,90,110,130,150,170,190,210,230,250,270,290,310,330)


yvalues <- c(msqcondfrail,msqcondgap,msqcondelapse,msqag,msqfrailelapse)
xvalues <- c(rep(sample,5))
names <- c(rep("Conditional Frailty",17), rep("Conditional, Gap",17), rep("Conditional, Elapse",17), rep("Andersen-Gill",17), rep("Frailty, Elapse",17))
ednohet <- data.frame(msqerror=yvalues,sample=xvalues,Model=names)

exp1 <- read.table(file="exp_noed_het_samplesize_ver210.txt", header=TRUE)
exp2 <- read.table(file="exp_noed_het_samplesize_ver230.txt", header=TRUE)
exp3 <- read.table(file="exp_noed_het_samplesize_ver250.txt", header=TRUE)
exp4 <- read.table(file="exp_noed_het_samplesize_ver270.txt", header=TRUE)
exp5 <- read.table(file="exp_noed_het_samplesize_ver290.txt", header=TRUE)
exp6 <- read.table(file="exp_noed_het_samplesize_ver2110.txt", header=TRUE)
exp7 <- read.table(file="exp_noed_het_samplesize_ver2130.txt", header=TRUE)
exp8 <- read.table(file="exp_noed_het_samplesize_ver2150.txt", header=TRUE)
exp9 <- read.table(file="exp_noed_het_samplesize_ver2170.txt", header=TRUE)
exp10 <- read.table(file="exp_noed_het_samplesize_ver2190.txt", header=TRUE)
exp11 <- read.table(file="exp_noed_het_samplesize_ver2210.txt", header=TRUE)
exp12 <- read.table(file="exp_noed_het_samplesize_ver2230.txt", header=TRUE)
exp13 <- read.table(file="exp_noed_het_samplesize_ver2250.txt", header=TRUE)
exp14 <- read.table(file="exp_noed_het_samplesize_ver2270.txt", header=TRUE)
exp15 <- read.table(file="exp_noed_het_samplesize_ver2290.txt", header=TRUE)
exp16 <- read.table(file="exp_noed_het_samplesize_ver2310.txt", header=TRUE)
exp17 <- read.table(file="exp_noed_het_samplesize_ver2330.txt", header=TRUE)


msqer11 <- mean((exp1$betafrailgap-(-1))^2, na.rm=TRUE)
msqer21 <- mean((exp2$betafrailgap-(-1))^2, na.rm=TRUE)
msqer31 <- mean((exp3$betafrailgap-(-1))^2, na.rm=TRUE)
msqer41 <- mean((exp4$betafrailgap-(-1))^2, na.rm=TRUE)
msqer51 <- mean((exp5$betafrailgap-(-1))^2, na.rm=TRUE)
msqer61 <- mean((exp6$betafrailgap-(-1))^2, na.rm=TRUE)
msqer71 <- mean((exp7$betafrailgap-(-1))^2, na.rm=TRUE)
msqer81 <- mean((exp8$betafrailgap-(-1))^2, na.rm=TRUE)
msqer91 <- mean((exp9$betafrailgap-(-1))^2, na.rm=TRUE)
msqer101 <- mean((exp10$betafrailgap-(-1))^2, na.rm=TRUE)
msqer111 <- mean((exp11$betafrailgap-(-1))^2, na.rm=TRUE)
msqer121 <- mean((exp12$betafrailgap-(-1))^2, na.rm=TRUE)
msqer131 <- mean((exp13$betafrailgap-(-1))^2, na.rm=TRUE)
msqer141 <- mean((exp14$betafrailgap-(-1))^2, na.rm=TRUE)
msqer151 <- mean((exp15$betafrailgap-(-1))^2, na.rm=TRUE)
msqer161 <- mean((exp16$betafrailgap-(-1))^2, na.rm=TRUE)
msqer171 <- mean((exp17$betafrailgap-(-1))^2, na.rm=TRUE)

msqer12 <- mean((exp1$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer22 <- mean((exp2$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer32 <- mean((exp3$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer42 <- mean((exp4$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer52 <- mean((exp5$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer62 <- mean((exp6$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer72 <- mean((exp7$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer82 <- mean((exp8$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer92 <- mean((exp9$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer102 <- mean((exp10$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer112 <- mean((exp11$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer122 <- mean((exp12$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer132 <- mean((exp13$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer142 <- mean((exp14$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer152 <- mean((exp15$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer162 <- mean((exp16$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer172 <- mean((exp17$betafrailelapse-(-1))^2, na.rm=TRUE)

msqer13 <- mean((exp1$betaag-(-1))^2, na.rm=TRUE)
msqer23 <- mean((exp2$betaag-(-1))^2, na.rm=TRUE)
msqer33 <- mean((exp3$betaag-(-1))^2, na.rm=TRUE)
msqer43 <- mean((exp4$betaag-(-1))^2, na.rm=TRUE)
msqer53 <- mean((exp5$betaag-(-1))^2, na.rm=TRUE)
msqer63 <- mean((exp6$betaag-(-1))^2, na.rm=TRUE)
msqer73 <- mean((exp7$betaag-(-1))^2, na.rm=TRUE)
msqer83 <- mean((exp8$betaag-(-1))^2, na.rm=TRUE)
msqer93 <- mean((exp9$betaag-(-1))^2, na.rm=TRUE)
msqer103 <- mean((exp10$betaag-(-1))^2, na.rm=TRUE)
msqer113 <- mean((exp11$betaag-(-1))^2, na.rm=TRUE)
msqer123 <- mean((exp12$betaag-(-1))^2, na.rm=TRUE)
msqer133 <- mean((exp13$betaag-(-1))^2, na.rm=TRUE)
msqer143 <- mean((exp14$betaag-(-1))^2, na.rm=TRUE)
msqer153 <- mean((exp15$betaag-(-1))^2, na.rm=TRUE)
msqer163 <- mean((exp16$betaag-(-1))^2, na.rm=TRUE)
msqer173 <- mean((exp17$betaag-(-1))^2, na.rm=TRUE)

msqer14 <- mean((exp1$betacondgap-(-1))^2, na.rm=TRUE)
msqer24 <- mean((exp2$betacondgap-(-1))^2, na.rm=TRUE)
msqer34 <- mean((exp3$betacondgap-(-1))^2, na.rm=TRUE)
msqer44 <- mean((exp4$betacondgap-(-1))^2, na.rm=TRUE)
msqer54 <- mean((exp5$betacondgap-(-1))^2, na.rm=TRUE)
msqer64 <- mean((exp6$betacondgap-(-1))^2, na.rm=TRUE)
msqer74 <- mean((exp7$betacondgap-(-1))^2, na.rm=TRUE)
msqer84 <- mean((exp8$betacondgap-(-1))^2, na.rm=TRUE)
msqer94 <- mean((exp9$betacondgap-(-1))^2, na.rm=TRUE)
msqer104 <- mean((exp10$betacondgap-(-1))^2, na.rm=TRUE)
msqer114 <- mean((exp11$betacondgap-(-1))^2, na.rm=TRUE)
msqer124 <- mean((exp12$betacondgap-(-1))^2, na.rm=TRUE)
msqer134 <- mean((exp13$betacondgap-(-1))^2, na.rm=TRUE)
msqer144 <- mean((exp14$betacondgap-(-1))^2, na.rm=TRUE)
msqer154 <- mean((exp15$betacondgap-(-1))^2, na.rm=TRUE)
msqer164 <- mean((exp16$betacondgap-(-1))^2, na.rm=TRUE)
msqer174 <- mean((exp17$betacondgap-(-1))^2, na.rm=TRUE)

msqer15 <- mean((exp1$betacondelapse-(-1))^2, na.rm=TRUE)
msqer25 <- mean((exp2$betacondelapse-(-1))^2, na.rm=TRUE)
msqer35 <- mean((exp3$betacondelapse-(-1))^2, na.rm=TRUE)
msqer45 <- mean((exp4$betacondelapse-(-1))^2, na.rm=TRUE)
msqer55 <- mean((exp5$betacondelapse-(-1))^2, na.rm=TRUE)
msqer65 <- mean((exp6$betacondelapse-(-1))^2, na.rm=TRUE)
msqer75 <- mean((exp7$betacondelapse-(-1))^2, na.rm=TRUE)
msqer85 <- mean((exp8$betacondelapse-(-1))^2, na.rm=TRUE)
msqer95 <- mean((exp9$betacondelapse-(-1))^2, na.rm=TRUE)
msqer105 <- mean((exp10$betacondelapse-(-1))^2, na.rm=TRUE)
msqer115 <- mean((exp11$betacondelapse-(-1))^2, na.rm=TRUE)
msqer125 <- mean((exp12$betacondelapse-(-1))^2, na.rm=TRUE)
msqer135 <- mean((exp13$betacondelapse-(-1))^2, na.rm=TRUE)
msqer145 <- mean((exp14$betacondelapse-(-1))^2, na.rm=TRUE)
msqer155 <- mean((exp15$betacondelapse-(-1))^2, na.rm=TRUE)
msqer165 <- mean((exp16$betacondelapse-(-1))^2, na.rm=TRUE)
msqer175 <- mean((exp17$betacondelapse-(-1))^2, na.rm=TRUE)

msqcondfrail <- c(msqer11, msqer21, msqer31, msqer41,msqer51,msqer61,msqer71, msqer81, msqer91, msqer101,msqer111, msqer121, msqer131, msqer141,msqer151,msqer161,msqer171)
msqag <-  c(msqer13, msqer23, msqer33, msqer43,msqer53,msqer63,msqer73, msqer83, msqer93,msqer103,msqer113, msqer123, msqer133, msqer143,msqer153,msqer163,msqer173)
msqfrailelapse <- c(msqer12, msqer22, msqer32, msqer42,msqer52,msqer62,msqer72, msqer82, msqer92,msqer102, msqer112, msqer122, msqer132, msqer142,msqer152,msqer162,msqer172)
msqcondgap <- c(msqer14, msqer24, msqer34, msqer44,msqer54,msqer64,msqer74, msqer84, msqer94, msqer104, msqer114, msqer124, msqer134, msqer144,msqer154,msqer164,msqer174)
msqcondelapse <- c(msqer15, msqer25, msqer35, msqer45,msqer55,msqer65,msqer75, msqer85, msqer95, msqer105, msqer115, msqer125, msqer135, msqer145,msqer155,msqer165,msqer175)
      
sample <- c(10,30,50,70,90,110,130,150,170,190,210,230,250,270,290,310,330)

yvalues <- c(msqcondfrail,msqcondgap,msqcondelapse,msqag,msqfrailelapse)
xvalues <- c(rep(sample,5))
names <- c(rep("Conditional Frailty",17), rep("Conditional, Gap",17), rep("Conditional, Elapse",17), rep("Andersen-Gill",17), rep("Frailty, Elapse",17))
noedhet <- data.frame(msqerror=yvalues,sample=xvalues,Model=names)

exp1 <- read.table(file="exp_ed_het_samplesize_ver210.txt", header=TRUE)
exp2 <- read.table(file="exp_ed_het_samplesize_ver230.txt", header=TRUE)
exp3 <- read.table(file="exp_ed_het_samplesize_ver250.txt", header=TRUE)
exp4 <- read.table(file="exp_ed_het_samplesize_ver270.txt", header=TRUE)
exp5 <- read.table(file="exp_ed_het_samplesize_ver290.txt", header=TRUE)
exp6 <- read.table(file="exp_ed_het_samplesize_ver2110.txt", header=TRUE)
exp7 <- read.table(file="exp_ed_het_samplesize_ver2130.txt", header=TRUE)
exp8 <- read.table(file="exp_ed_het_samplesize_ver2150.txt", header=TRUE)
exp9 <- read.table(file="exp_ed_het_samplesize_ver2170.txt", header=TRUE)
exp10 <- read.table(file="exp_ed_het_samplesize_ver2190.txt", header=TRUE)
exp11 <- read.table(file="exp_ed_het_samplesize_ver2210.txt", header=TRUE)
exp12 <- read.table(file="exp_ed_het_samplesize_ver2230.txt", header=TRUE)
exp13 <- read.table(file="exp_ed_het_samplesize_ver2250.txt", header=TRUE)
exp14 <- read.table(file="exp_ed_het_samplesize_ver2270.txt", header=TRUE)
exp15 <- read.table(file="exp_ed_het_samplesize_ver2290.txt", header=TRUE)
exp16 <- read.table(file="exp_ed_het_samplesize_ver2310.txt", header=TRUE)
exp17 <- read.table(file="exp_ed_het_samplesize_ver2330.txt", header=TRUE)


msqer11 <- mean((exp1$betafrailgap-(-1))^2, na.rm=TRUE)
msqer21 <- mean((exp2$betafrailgap-(-1))^2, na.rm=TRUE)
msqer31 <- mean((exp3$betafrailgap-(-1))^2, na.rm=TRUE)
msqer41 <- mean((exp4$betafrailgap-(-1))^2, na.rm=TRUE)
msqer51 <- mean((exp5$betafrailgap-(-1))^2, na.rm=TRUE)
msqer61 <- mean((exp6$betafrailgap-(-1))^2, na.rm=TRUE)
msqer71 <- mean((exp7$betafrailgap-(-1))^2, na.rm=TRUE)
msqer81 <- mean((exp8$betafrailgap-(-1))^2, na.rm=TRUE)
msqer91 <- mean((exp9$betafrailgap-(-1))^2, na.rm=TRUE)
msqer101 <- mean((exp10$betafrailgap-(-1))^2, na.rm=TRUE)
msqer111 <- mean((exp11$betafrailgap-(-1))^2, na.rm=TRUE)
msqer121 <- mean((exp12$betafrailgap-(-1))^2, na.rm=TRUE)
msqer131 <- mean((exp13$betafrailgap-(-1))^2, na.rm=TRUE)
msqer141 <- mean((exp14$betafrailgap-(-1))^2, na.rm=TRUE)
msqer151 <- mean((exp15$betafrailgap-(-1))^2, na.rm=TRUE)
msqer161 <- mean((exp16$betafrailgap-(-1))^2, na.rm=TRUE)
msqer171 <- mean((exp17$betafrailgap-(-1))^2, na.rm=TRUE)

msqer12 <- mean((exp1$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer22 <- mean((exp2$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer32 <- mean((exp3$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer42 <- mean((exp4$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer52 <- mean((exp5$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer62 <- mean((exp6$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer72 <- mean((exp7$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer82 <- mean((exp8$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer92 <- mean((exp9$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer102 <- mean((exp10$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer112 <- mean((exp11$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer122 <- mean((exp12$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer132 <- mean((exp13$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer142 <- mean((exp14$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer152 <- mean((exp15$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer162 <- mean((exp16$betafrailelapse-(-1))^2, na.rm=TRUE)
msqer172 <- mean((exp17$betafrailelapse-(-1))^2, na.rm=TRUE)

msqer13 <- mean((exp1$betaag-(-1))^2, na.rm=TRUE)
msqer23 <- mean((exp2$betaag-(-1))^2, na.rm=TRUE)
msqer33 <- mean((exp3$betaag-(-1))^2, na.rm=TRUE)
msqer43 <- mean((exp4$betaag-(-1))^2, na.rm=TRUE)
msqer53 <- mean((exp5$betaag-(-1))^2, na.rm=TRUE)
msqer63 <- mean((exp6$betaag-(-1))^2, na.rm=TRUE)
msqer73 <- mean((exp7$betaag-(-1))^2, na.rm=TRUE)
msqer83 <- mean((exp8$betaag-(-1))^2, na.rm=TRUE)
msqer93 <- mean((exp9$betaag-(-1))^2, na.rm=TRUE)
msqer103 <- mean((exp10$betaag-(-1))^2, na.rm=TRUE)
msqer113 <- mean((exp11$betaag-(-1))^2, na.rm=TRUE)
msqer123 <- mean((exp12$betaag-(-1))^2, na.rm=TRUE)
msqer133 <- mean((exp13$betaag-(-1))^2, na.rm=TRUE)
msqer143 <- mean((exp14$betaag-(-1))^2, na.rm=TRUE)
msqer153 <- mean((exp15$betaag-(-1))^2, na.rm=TRUE)
msqer163 <- mean((exp16$betaag-(-1))^2, na.rm=TRUE)
msqer173 <- mean((exp17$betaag-(-1))^2, na.rm=TRUE)

msqer14 <- mean((exp1$betacondgap-(-1))^2, na.rm=TRUE)
msqer24 <- mean((exp2$betacondgap-(-1))^2, na.rm=TRUE)
msqer34 <- mean((exp3$betacondgap-(-1))^2, na.rm=TRUE)
msqer44 <- mean((exp4$betacondgap-(-1))^2, na.rm=TRUE)
msqer54 <- mean((exp5$betacondgap-(-1))^2, na.rm=TRUE)
msqer64 <- mean((exp6$betacondgap-(-1))^2, na.rm=TRUE)
msqer74 <- mean((exp7$betacondgap-(-1))^2, na.rm=TRUE)
msqer84 <- mean((exp8$betacondgap-(-1))^2, na.rm=TRUE)
msqer94 <- mean((exp9$betacondgap-(-1))^2, na.rm=TRUE)
msqer104 <- mean((exp10$betacondgap-(-1))^2, na.rm=TRUE)
msqer114 <- mean((exp11$betacondgap-(-1))^2, na.rm=TRUE)
msqer124 <- mean((exp12$betacondgap-(-1))^2, na.rm=TRUE)
msqer134 <- mean((exp13$betacondgap-(-1))^2, na.rm=TRUE)
msqer144 <- mean((exp14$betacondgap-(-1))^2, na.rm=TRUE)
msqer154 <- mean((exp15$betacondgap-(-1))^2, na.rm=TRUE)
msqer164 <- mean((exp16$betacondgap-(-1))^2, na.rm=TRUE)
msqer174 <- mean((exp17$betacondgap-(-1))^2, na.rm=TRUE)

msqer15 <- mean((exp1$betacondelapse-(-1))^2, na.rm=TRUE)
msqer25 <- mean((exp2$betacondelapse-(-1))^2, na.rm=TRUE)
msqer35 <- mean((exp3$betacondelapse-(-1))^2, na.rm=TRUE)
msqer45 <- mean((exp4$betacondelapse-(-1))^2, na.rm=TRUE)
msqer55 <- mean((exp5$betacondelapse-(-1))^2, na.rm=TRUE)
msqer65 <- mean((exp6$betacondelapse-(-1))^2, na.rm=TRUE)
msqer75 <- mean((exp7$betacondelapse-(-1))^2, na.rm=TRUE)
msqer85 <- mean((exp8$betacondelapse-(-1))^2, na.rm=TRUE)
msqer95 <- mean((exp9$betacondelapse-(-1))^2, na.rm=TRUE)
msqer105 <- mean((exp10$betacondelapse-(-1))^2, na.rm=TRUE)
msqer115 <- mean((exp11$betacondelapse-(-1))^2, na.rm=TRUE)
msqer125 <- mean((exp12$betacondelapse-(-1))^2, na.rm=TRUE)
msqer135 <- mean((exp13$betacondelapse-(-1))^2, na.rm=TRUE)
msqer145 <- mean((exp14$betacondelapse-(-1))^2, na.rm=TRUE)
msqer155 <- mean((exp15$betacondelapse-(-1))^2, na.rm=TRUE)
msqer165 <- mean((exp16$betacondelapse-(-1))^2, na.rm=TRUE)
msqer175 <- mean((exp17$betacondelapse-(-1))^2, na.rm=TRUE)

msqcondfrail <- c(msqer11, msqer21, msqer31, msqer41,msqer51,msqer61,msqer71, msqer81, msqer91, msqer101,msqer111, msqer121, msqer131, msqer141,msqer151,msqer161,msqer171)
msqag <-  c(msqer13, msqer23, msqer33, msqer43,msqer53,msqer63,msqer73, msqer83, msqer93,msqer103,msqer113, msqer123, msqer133, msqer143,msqer153,msqer163,msqer173)
msqfrailelapse <- c(msqer12, msqer22, msqer32, msqer42,msqer52,msqer62,msqer72, msqer82, msqer92,msqer102, msqer112, msqer122, msqer132, msqer142,msqer152,msqer162,msqer172)
msqcondgap <- c(msqer14, msqer24, msqer34, msqer44,msqer54,msqer64,msqer74, msqer84, msqer94, msqer104, msqer114, msqer124, msqer134, msqer144,msqer154,msqer164,msqer174)
msqcondelapse <- c(msqer15, msqer25, msqer35, msqer45,msqer55,msqer65,msqer75, msqer85, msqer95, msqer105, msqer115, msqer125, msqer135, msqer145,msqer155,msqer165,msqer175)
      
sample <- c(10,30,50,70,90,110,130,150,170,190,210,230,250,270,290,310,330)

yvalues <- c(msqcondfrail,msqcondgap,msqcondelapse,msqag,msqfrailelapse)
xvalues <- c(rep(sample,5))
names <- c(rep("Conditional Frailty",17), rep("Conditional, Gap",17), rep("Conditional, Elapse",17), rep("Andersen-Gill",17), rep("Frailty, Elapse",17))
edhet <- data.frame(msqerror=yvalues,sample=xvalues,Model=names)

msqplot1 <- qplot(sample,msqerror,data=noednohet,geom="line",colour=Model,linetype=Model) + scale_colour_brewer(pal="Dark2") + labs(y="Mean Squared Error", x="Cases (Maximum Events=3)") + opts(legend.position="none") + scale_linetype_manual("Model",values=c(1:2,4:6)) + opts(title="No Event Dependence, No Heterogeneity") +ylim(0,.6)

msqplot2 <- qplot(sample,msqerror,data=ednohet,geom="line",colour=Model,linetype=Model) + scale_colour_brewer(pal="Dark2") + labs(y="Mean Squared Error", x="Cases (Maximum Events=3)", main="Event Dependence, No Heterogeneity") + opts(legend.position=c(.7,.7)) + scale_linetype_manual("Model",values=c(1:2,4:6))  + opts(title="Event Dependence, No Heterogeneity") 

msqplot3 <- qplot(sample,msqerror,data=noedhet,geom="line",colour=Model,linetype=Model) + scale_colour_brewer(pal="Dark2") + labs(y="Mean Squared Error", x="Cases (Maximum Events=3)", main="No Event Dependence, Heterogeneity") + opts(legend.position="none")  + scale_linetype_manual("Model",values=c(1:2,4:6))  + opts(title="No Event Dependence, Heterogeneity") 

msqplot4 <- qplot(sample,msqerror,data=edhet,geom="line",colour=Model,linetype=Model) + scale_colour_brewer(pal="Dark2") + labs(y="Mean Squared Error", x="Cases (Maximum Events=3)", main="Event Dependence, Heterogeneity") + opts(legend.position="none")  + scale_linetype_manual("Model",values=c(1:2,4:6))  + opts(title="Event Dependence, Heterogeneity")


## B & W


msqplot1 <- qplot(sample,msqerror,data=noednohet,geom="line",linetype=Model) + theme_bw()+ labs(y="Mean Squared Error", x="Cases (Maximum Events=3)") + opts(legend.position="none") + scale_linetype_manual("Model",values=c(1,2,4,3,5)) + opts(title="No Event Dependence, No Heterogeneity") +ylim(0,.6)

msqplot2 <- qplot(sample,msqerror,data=ednohet,geom="line",linetype=Model) + theme_bw() + labs(y="Mean Squared Error", x="Cases (Maximum Events=3)", main="Event Dependence, No Heterogeneity") + opts(legend.position=c(.7,.7)) + scale_linetype_manual("Model",values=c(1,2,4,3,5))  + opts(title="Event Dependence, No Heterogeneity") 

msqplot3 <- qplot(sample,msqerror,data=noedhet,geom="line",linetype=Model) + theme_bw() + labs(y="Mean Squared Error", x="Cases (Maximum Events=3)", main="No Event Dependence, Heterogeneity") + opts(legend.position="none")  + scale_linetype_manual("Model",values=c(1,2,4,3,5))  + opts(title="No Event Dependence, Heterogeneity") 

msqplot4 <- qplot(sample,msqerror,data=edhet,geom="line",linetype=Model) + theme_bw()+ labs(y="Mean Squared Error", x="Cases (Maximum Events=3)", main="Event Dependence, Heterogeneity") + opts(legend.position="none")  + scale_linetype_manual("Model",values=c(1,2,4,3,5))  + opts(title="Event Dependence, Heterogeneity")


pdf("mse_maxn_panel.pdf",height=8.5,width=11)
grid.newpage()
pushViewport(viewport(layout = grid.layout(2,2)))
vplayout <- function(x,y)
viewport(layout.pos.row=x,layout.pos.col=y)
print(msqplot1,vp=vplayout(1,1))
print(msqplot2,vp=vplayout(1,2))
print(msqplot3,vp=vplayout(2,1))
print(msqplot4,vp=vplayout(2,2))
dev.off()
